////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -------------------- Ebrahim Foulaadvand & Somayyeh Belbasi, 23 June 2009 --------------------------- //
//                                                                                                        //
//  The programme "CoupledOscillators" evaluates the temporal evolution of a system of coupled harmonic   //
//  oscillator which are connected to each other by springs. The fix boundary condition is used.          //
//  The second order Runge-Kutta algorithm is used for solving the motion equations.                      //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

main()

{

randomize();

const N=10, T=5000; // N= number of oscillators, T= number of simulation timesteps.

double delt=0.01; // delt=timestep.

int i,j,n,s,t;

vector< double > m(N+2,0),k(N+3,0),xanal(N+2,0),vanal(N+2,0),x(N+2,0),v(N+2,0),k1x(N+2,0),k2x(N+2,0),
                 k1v(N+2,0),k2v(N+2,0),nu(N+2,0),mu(N+2,0),w(N+2,0);

// Arrays "m" and "k" store the oscillators masses and the spring constants. Arrays "x" and "v" store
// oscillators positions and velocities. Arrays "xanal" and "vanal" store the analytical values of
// positions and velocitis. Arrays "nu","mu" and "w" are used for analytical solutions. Refer to book text
// to see their definitions.


ofstream file1 ("x analytic5.plt"); // output file for the analytic poistion of particle 5.
ofstream file2 ("v analytic5.plt"); // output file for the analytic velocity of particle 5.
ofstream file3 ("phasexv analytic5.plt"); // output file for the analytic phase space plot of particle 5.

ofstream file4 ("x5.plt"); // output file for the computed poistion of particle 5.
ofstream file5 ("v5.plt"); // output file for the computed velocity of particle 5.
ofstream file6 ("phasexv5.plt"); // output file for the computed phase space plot of particle 5.


for (int i=1;i<=N;i++){ m[i]=1; k[i]=1; } // specification of masses and spring values.

k[N+1]=1;


// ---------------------------------------- Initial Conditions  --------------------------------------------------

xanal[0]=0; vanal[0]=0;
xanal[N+1]=0; vanal[N+1]=0;

//The fixed boundary condition is modeled by introducing additional fixed particles 0 and N+1.

for (int i=1;i<=N;i++){ xanal[i]=0; vanal[i]=0; }

x[1]=1;

xanal[1]=1;


//----------------------------------------- Analytical Solution ------------------------------------------


for (s=1;s<=N;s++) { mu[s]=0; nu[s]=0; w[s]=2*sqrt(k[s]/m[s])*sin( 0.5*s*M_PI/(N+1) ); }

for (s=1;s<=N;s++) {

for (j=1;j<=N;j++) { mu[s]+=(2/double (N+1))*xanal[j]*sin(j*M_PI*s/(N+1)); nu[s]+=(-2/(w[s]*(N+1)))*vanal[j]*sin(j*M_PI*s/(N+1)); }

//cout<<"mu["<<s<<"]= "<<mu[s]<<endl;//getch();

//cout<<"nu["<<s<<"]= "<<nu[s]<<endl;getch();

}

for (int t=1;t<=T;t++) { //cout<<"t= "<<t<<endl;getch();


for (j=1;j<=N;j++) { xanal[j]=0; vanal[j]=0;

for (n=1;n<=N;n++) {

xanal[j]+=sin(j*M_PI*n/double (N+1))*( mu[n]*cos(w[n]*t*delt) - nu[n]*sin(w[n]*t*delt));

vanal[j]+=sin(j*M_PI*n/double (N+1))*( -mu[n]*w[n]*sin(w[n]*t*delt) - nu[n]*w[n]*cos(w[n]*t*delt));

}

} // j

//cout<<"x["<<5<<"]= "<<xanal[5]<<endl;getch();

file1<<t*delt<<" "<<xanal[5]<<endl; // storing in the data file
file2<<t*delt<<" "<<vanal[5]<<endl; // storing in the data file
file3<<xanal[5]<<" "<<vanal[5]<<endl; // storing in the data file


} // t


//-----------  2nd Order Runge-Kutta algorithm --------------------

x[0]=0; v[0]=0;

x[N+1]=0; v[N+1]=0;

for (int i=1;i<=N;i++){x[i]=0; v[i]=0; }

x[1]=1.0;

for (int t=1;t<=T;t++) { //cout<<"t= "<<t<<endl;getch();


for (i=1;i<=N;i++) {

k1v[i]=-delt*(k[i]/m[i])*(x[i]-x[i-1])-delt*(k[i+1]/m[i])*(x[i]-x[i+1]);
k1x[i]=delt*v[i];


k2v[i]=-delt*(k[i]/m[i])*(x[i]+0.5*k1x[i]-x[i-1]-0.5*k1x[i-1])-delt*(k[i+1]/m[i])*(x[i]+0.5*k1x[i]-x[i+1]-0.5*k1x[i+1]);
k2x[i]=delt*(v[i] +0.5*k1v[i]);


v[i]=v[i] + k2v[i];
x[i]=x[i] + k2x[i];

} // i

//cout<<"x["<<9<<"]= "<<x[9]<<endl;getch();

file4<<t*delt<<" "<<x[5]<<endl; // storing in the data file
file5<<t*delt<<" "<<v[5]<<endl; // storing in the data file
file6<<x[5]<<" "<<v[5]<<endl;   // storing in the data file


} // t


cout<<"end."<<endl;getch();


return 0;

} //main